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Abstract 

This  paper  is  concerned  with  the  identification  of  the  geometrical  structure  of  the 
system  boundary  for  a  two-dimensional  diffusion  system.  The  domain  identification 
problem  treated  here  is  converted  into  an  optimization  problem  based  on  a  fit-to-data 
criterion  and  theoretical  convergence  results  for  approximate  identification  techniques 
are  discussed.  Results  of  numerical  experiments  to  demonstrate  the  efficacy  of  the 
theoretical  ideas  are  reported. 
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I.  INTRODUCTION 


'  Domain  identification  problems  are  important  in  the  design  of  engineering  systems  and 
frequently  such  problems  are  treated  as  a  branch  of  the  calculus  of  variations  which  involves 
nonlinear  optimization  techniques,  optimal  control  theory,  partial  differential  equations  (el¬ 
liptic,  parabolic,  hyperbolic,  etc.)  and  related  numerical  methods.  Domain  identification 
for  elliptic  systems  has  been  studied  theoretically  and  numerically  by  many  authors  (see 
e.g.,  [5],[7],[lO],[l3]).  For  parabolic  systems,  a  couple  of  numerical  methods  for  identi¬ 
fying  the  domain  or  boundary  have  been  investigated. in  [14], [15).  -  Until  recently,  most 
investigations  concentrated  on  the  “optimal  shape  design  problem”  which  is  motivated  by 

numerous  applications  to  structural,  engine,  airplane,  ship  designs,  etc.  (see  [10]  and  the 

\ 

references  therein). ^n  this  paper,  our  concern  for  domain  identification  is  motivated  by  an 

application  that  is  different  from  these  shape  design  problems.  However,  as  we  shall  see, 

the  resulting  theoretical  aspects  are  closely  related.  Recently,  associated  with  the  use  of 

fiber  reinforced  composite  materials  for  aerospace  structures,  there  is  growing  interest  in 

the  detection  and  characterization  of  large  structural  flaws  which  may  not  be  detectable 

by  visual  inspection.  One  recent  effort  has  focused  on  non-destructive  evaluation  methods 

(NDE)  based  on  the  measurement  of  thermal  diffusivity  in  composite  materials  (see  e.g., 

[8]).  Motivated  by  these  problems,  we  consider  the  domain  identification  for  parabolic 

systems.  - 
l 

To  explain  our  approach,  we  restrict  our  attention  to  a  2-D  domain  identification  prob¬ 
lem.  We  consider  the  bounded  domain  G(q)  in  two-dimensional  Euclidean  space  as  follows: 

G(q)  =  {(x!,2:2)|0  <  xi  <  1,  0  <  x2  <  r(xi,<?)} 

where  Xi  — »  r(ii,q)  is  some  parameterized  real  function  which  is  assumed  to  characterize 
the  unknown  part  of  the  boundary  and  q  is  a  constant  parameterization  vector  to  be 
identified  among  values  in  a  given  compact  admissible  parameter  set  Q.  As  depicted  in 


Fig.  1.1.  The  spatial  domain  G  and  its  boundary  dG j,  dG j,  dGq,  dGt. 


Fig.  1.1,  we  assume  the  boundary  of  G(g)  consists  of  the  following  components: 
dG i  —  {x  =  (zi,i2)|0  <  X\  <  1,  x2  =  0} 

dG2  =  {x  =  ^Xi ,  X2)  jxi  =  1,  0  <  x2  <  t} 

dGq  =  {x  =  (x!,X2)|0  <  xi  <  1,  x2  =  r(x!,9)} 

dGi  —  {x  =  (x1,X2)jxi  =  0,  0  <  x2  <  1} 

The  measurement  system  is  described  by  the  following  2-D  diffusion  equation: 

—  1 1  ’  —  —  CiAu(t,x)  +  c0u(t,x)  =  0  in  T  x  G(q)  ( 

at 

with  the  initial  and  boundary  conditions, 


tt(0,  x)  =  £o(x)  on  G{q) 
2 


c  v  r 


KJJX 

—ci— - h  h(u  —  /)  =  0  on  T  x  3Gi 

OX  2 


du 

ci  — - (-  hu  =  0 

uXi 

du 
dn  = 
du, 

-Cj- - 1-  Au  =  0 


on  T  x  dG2 


on  T  x  dG0 


on  T  x  <9G4 


(I.lc) 

(l.W) 

(lit) 

(11/) 


where  Cj,  Co  and  A  are  thermal  diffusivity,  radiation  coefficient  and  heat  transfer  coefficient, 
respectively,  which  are  given  constants,  and  where  T  denotes  the  time  interval  (0,  tf)  during 
which  the  process  is  observed.  In  the  above  system,  /  is  the  known  boundary  input  defined 
on  T  x  dG\  and  uq  is  the  given  initial  function  defined  on  Q  where  O  is  a  known  bounded 
domain  in  R 2  such  that  f!  D  G(q)  for  any  qeQ.  The  system  output  is  assumed  to  be  on  a 
subset  £  of  the  boundary  dG\,  and  mathematically,  the  observation  is  taken  as 

y(t,x,q)  =  u{t,xu0,q)  (t,x)eT  x  £.  (1.2) 

From  a  physical  point  of  view,  the  system  state  u  =  u(f,x)  represents  the  temperature 
distribution  at  time  t  at  location  x  =  (x^  x2)  and,  the  boundary  input  /  and  the  output 
y  correspond,  respectively,  to  the  thermal  source  (for  example,  by  a  laser  beam)  and 
the  observation  of  the  temperature  distribution  at  the  surface  of  the  material  (e.g.,  by  an 
infrared  imager)  (see  [8]  for  more  details).  Thus,  the  search  for  structural  flaws  in  materials 
may  be  formulated  as  an  inverse  problem  for  a  heat  diffusion  system.  The  problem  treated 
here  is  that  of  identifying,  from  input  and  output  data  {/,  u0,y}  on  ( T  x  dG j)  x  G  x  (T  x 
£),  the  constant  parameter  vector  q  in  Q  determining  the  geometrical  structure  of  the 
boundary  dGq. 

In  Section  2,  we  formulate  this  problem  in  an  abstract  setting  in  a  Hilbert  space.  In  Sec¬ 
tion  3,  for  computational  purposes,  we  approximate  the  Hilbert  space  by  finite  dimensional 
subspaces  and  we  discuss  the  convergence  analysis  for  the  approximate  identification  prob¬ 
lems.  In  Section  4,  a  practical  optimization  technique  based  on  a  finite  element  approach 
is  outlined.  Some  numerical  results  for  a  simple  example  are  given  in  Section  5. 
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II.  PROBLEM  FORMULATION  AND  BASIC  ASSUMPTIONS 


For  the  discussions  here,  we  restrict  the  geometrical  structure  of  the  boundary  dGq  by 
imposing  the  following  hypotheses: 

(H-0)  The  admissible  parameter  set  Q  is  a  compact  subset  of  Rn; 

(H-l)  For  each  qeQ,  we  have  r(<7)eW^(0, 1); 

(H-2)  For  each  qeQ,  we  have  r(0,  q)  —  r(l,  q)  =  t; 

(H-3)  There  are  constants  Pi  and  P2  satisfying  0<  Pi  <  i  <  P2  <  oo  such  that,  for  qeQ, 
we  have 

Pi<r{£,q)<P2  a.e.  in  (0,1); 

and 

(H-4)  There  exists  a  constant  M  such  that 

K£,9)  -r(£,g)|1>00  <  M\q-q\  for  q,  qeQ,  £e[0,  l], 

where  |  ■  |i(00  denotes  the  norm  of  VV^O,  1). 

We  make  the  following  assumptions  for  the  class  of  system  inputs: 

(H-5)  u0eL2(  fi); 

and 

(H-6)  feL2(T ;  L2(0, 1)). 

It  follows  from  results  in  ([9],  Ch.  3)  that,  under  the  hypotheses  (H-l)-(H-3),  (H-5),  and 
(H-6),  for  each  fixed  qeQ,  there  exists  for  (1.1)  a  unique  solution  u  in  L2{T\  H1{G(q))). 
Following  standard  procedures  in  optimal  shape  design  techniques  ([10],  Ch.  8,  p.  125), 
we  introduce  the  affine  mapping 


(zi,e2)  =  T{q){xi,x2) 


given  by 


z\  =  zi 

_  tx-i 
r{xuq)' 

Note  that  this  is  equivalent  to  x\  =  Zi,x2  =  r(zi,q)z2/ £.  Under  this  coordinate  change,  the 
system  domain  G(q)  is  transformed  into  the  fixed  domain  G 

G{q)  — *G  =  (0,1)  x  (0,  £) 

which  is  independent  of  the  parameter  q.  Using  this  coordinate  transformation,  we  obtain 
the  system  state  u  given  by 

u(t,z)  =  u(t)  o  T~1(q)  =  u^x^Z!),  x2 (21,22)); 


this  transformed  state  then  satisfies  the  system  equation 


with 


u(0,  z)  =  u0(z)  o  T^1(q)  on  G 


in  T  x  G 
(2.1  a) 

(2.16) 


Cit  du 

r{z\,q)  dz2 


+  h{u  —  /)  =  0  on  T  x  dGi 


(2.1c) 


Cl'P~~T z1r'{l,q)^L  +  hu  =  0  on  T  x  dGi 

OZ\  t  C/Z2 


(2. Id) 


+  =  °  “  Txda 3  (21e) 


clp-  +  C-±z2r'(0,q)p^  +  hu  =  0  on  T  x  dGV 

c/^i  £  c/^2 


(2.1/) 


In  this  system  the  coefficients  are  given  by 


au  :=  ci 

/  \  /  \  c1r'(zuq)z2 

«»W  =  »,.(«)  !=  -  r(-t) 


(2'5) 

**w  =  (2-6) 
where  r'  denotes  dr jdz\,  while  dGq  has  been  mapped  into 

dGz  =  {z  —  (2l,22)|0  <  z\  <  1,  22  =  0- 

If  we  consider  a  variational  formulation  similar  to  that  in  [2]  ,[9],  the  system  dynamics  can 
be  described  by  the  variational  form: 

<  ~dP~’<t>  >  +cr(?M“(*M)  =  L{t,q){<t>)  for  4>eHl{G) 


u(0)  =  u0  o  T  *(g) 

where  the  bracket  <  V  >  denotes  the  scalar  product  in  L2(G)  and  where  cr(^)(-,-)  and 
L(q)(-)  denote,  respectively,  a  sesquilinear  form  on  Hl{G )  x  Hl{G )  and  a  linear  functional 
on  Hl(G).  Explicitly,  cr(g)(-,-)  and  L(g)(-)  are  given  by  for  <j>,xl)tHl{G)  by 


o(q)M)  :=  //.[  ^a,,(2)  — —  +  +  co U}dz  +  l  [U}dGl  dz> 

(2.8) 

+  JQ  hlW}dG3UdGtdz2’ 

f1  lh 

L{t,q)(<t>):=  - - rf{t,Zi)[<f>}dG,dzu  (2.9) 


amj$jKpy’m*h‘ 

respectively.  We  note  that,  from  (H-0)  to  (H-4), 


sup  \r'(z,q)\  <  R 

where  R  is  some  constant  independent  of  q.  Hence,  we  obtain 


\bj{z,q)\  <  K2  <  oo  for  j-  1,2, 


where  (assuming  R  >  1) 


K ,  = 


c^R2 


(2.18) 


(2.19) 


For  the  last  two  boundary  integrals  in  (2.8),  by  virtue  of  (H-3),  the  following  inequality 
holds: 

r 1  rt  r 

Jo  r(zi  ~q)  ^aG'dzi  +  J0  >  h  J^[cf>2}a^ds  '2.20) 

where  ds  denotes  a  line  element  on  dG  and  dG  —  dGi  U  dG2  U  dG4.  From  (2.14),  (2.18), 
and  (2.20),  we  can  derive  the  coercivity  property  of  the  sesquilinear  form.  Namely,  the 
sesquilinear  form  satisfies 

'<*><*•  « *  *  /  L  g 1 - *  /  L  g 

+h 

+  A’!  (//dgl3^|!‘i2  +  /ao^!|a5‘ls)  (2'21) 


where 


Ks  =  min  ■ 

Friedrich’s  second  inequality  ([l],  p.  124)  asserts  that  if  dG  is  a  nontrivial  subset  of  dG, 
then  there  exists  a  positive  constant  a  such  that 


/  4  g  ^ 1  Hz + U^ds  a  a{d) 


(2.22) 


fci: 


i 


K 
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By  applying  this  to  the  last  parenthesis  of  (2.21),  we  can  conclude  that 


where 


respectively. 


o(q)(<f>,<f>)  >  K4\<f>\2v  - 


To  prove  the  boundedness  of  a(q ),  we  note  that 


+  1  YlJ  i  +  l  / 

r1  ih  re 

+  1  /  "7 - r[<A0]aGj^i  +  / 

-'0  r^,?)  /o 


The  first  three  integrals  of  RHS  in  (2.24)  satisfy 

,  v-  f  f  .  /_  ^  ^  , 


^  f  faaii(Z'q}jz-dz-dz\-4  fU£  Wvl^lv 


*€{0,i]x[0,£] 


l5l/  /-  <  2  sup  (&>(*>  ?)  I  |0|v#|v  (2-: 

|  J  J  CQ<f>r}jdz\  <  c0\4>\v\4>\v ,  (2: 

respectively.  From  (2.2)  to  (2.4),  it  follows  that,  under  hypotheses  (H-l)  and  (H-3)  ( 
(2-17)), 


where 


sup  |a,y(z,g)|  <  Ke,  <  oo 

i,;<2 

Jf<(0, 1  j  X  (0,£j 


K6=^(fl’+1). 


From  (2.18)  and  (2.25)-(2.28),  we  have 


yi  ih  r( 

T{q){4>, VO!  £  (4*6  +  2/f2  -f  co)Hk|0|k  +  I  I  ~7 — — rfVV’JaGidzi  +  J  h[<f>ip\aGjUSG4dz 


Furthermore,  the  boundary  integral  term  satisfies 


r  1  £h  ri 

1/  -7 - A<t>iJj}dGidz1+  h{4>ip}aG2udG<dz2\  <  K7\<j>\v\4)\v 

Jo  r(zi,q)  Jo 

which  follows  from  the  trace  inequality  ( [l J ,  p.  124) 

fj<t>2\dGds<l(GM2v, 

J  oG 


where 


K7  = 


M^[G) 

Wi 


Consequently,  we  can  prove  the  boundedness  property  of  c(<7)(-,  •). 

To  establish  the  continuity  property,  we  note  that,  for  any  q  and  qeQ, 

-  o{q){<p,tp)\ 

+  1  /  /- 

J  ->G  j  OZj 

+l/o‘tt(i4^4))  lWla°1<i211 

+  J  fc  " MSI  lf~l  +  f  fd  IMS  -  MSI  l|j)l  Md; 

Under  the  hypotheses  (H-l)  and  (H-3),  we  argue  that 

c  £ 

|au(?)  ~  012(9)!  <  ^maz(^’  1)lr(9)  -  r(9) !j,oo 

,  ,  .  ....^2cxeR  ,2Rfh 

1 022(9)  -  022(9)1  <  — -M—‘‘ maX[—M-’  1)lr(9)  -  r (9)  11,00 


01 


PI 


|t»i(g)  -  61(9)1  <  J-^ax{^~,\)\r{q)  -  r(g)|liD 
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and 


M$)  -Mg)l  <  -^-rnax{-^-,l)\r{q)  -  r{q) j,,0 


l^y-4j!sil|r(?)-r®1|o(°'I»- 


Applying  these  inequalities  into  (2.32),  we  have 


k(g)(<£>0)  ~  *(§)(&  0)1  <  Kt\ r(g)  -  r(5)|itOo|^|v|0|v 


(2.33) 


where 


is  2  Cl£  ,J2  ^  2cii2  R  ,2Rfa  ^  ,  c,  ,  i2  ^,2  Cl£iZ  ,iZft  ,  th 


From  the  hypotheses  (H-4),  we  can  thus  infer  the  continuity  of  the  sesquilinear  form 
cr(g)(-,  •)  with  respect  to  the  parameter  q  in  Q.  The  proof  has  been  completed. 

For  the  system  (2.7),  the  output  can  be  represented  as  the  restriction  of  u[t)  to  a  subset 


E  C  dGx  of  positive  measure,  i.e, 


y{t,q)  =  u(*,g)|s. 


(2.34) 


We  assume  (see  (H-0))  throughout  that  the  admissible  parameter  set  Q  is  a  given  compact 


subset  of  Rn.  The  fundamental  identification  problem  considered  here  is  based  on  the 


fit-to-data  functional  (see  [2])  given  by 


J(g)  =  \  //  l|y(*>g)  -  MOIfc* 


(2.35) 


where  F  =  Z/2(E),  {yd{t)}uT  are  given  observed  data,  and  y(t,q )  is  the  solution  of  (2.7) 


corresponding  to  qeQ,  Then  our  problem  is  stated  as  follows: 


(IDP)  Find  q'eQ  which  minimizes  J(q)  given  in  (2.35)  subject  to  the  system  (2.7)  and 
(2.34). 


In  the  next  section,  we  consider  a  family  of  approximating  identification  problems  associ¬ 


ated  with  (IDP). 


III.  APPROXIMATE  IDENTIFICATION  PROBLEMS 


The  approximation  scheme  we  have  employed  is  based  on  the  use  of  a  finite  element 
Galerkin  approach  to  construct  a  sequence  of  finite  dimensional  approximating  identifica¬ 
tion  problems.  Let  us  choose  Ujv=i{$fr}£[:i  as  a  set  °f  basis  functions  in  Hl(G).  That  is, 
for  all  N,  are  linearly  independent  and  Uw  spanf^}^  is  dense  in  the  V  norm 

in  V  =  Hl{G).  We  choose  the  approximation  subspaces  as 


Hn  :=  span{(f>i  <£$}• 

Then,  we  can  define  the  approximate  solution  of  Eq.  (2.7)  by 


uN[t,q)  =  wif(t,‘l)<t>? 

i<N 

(3.1) 

where  u>^[t,  q) 

are  chosen  such  that  for  j  =  1, 2,  •  •  • ,  TV, 

<  dif,£q) ,*?  >  +*(»)(«"(<,«),*?)  - 

(3.2a) 

and 

«*(o)  =  <  u° 0 

i<N 

(3.26) 

Hence  the  system  (2.7)  and  the  output  (2.34)  can  be  approximated  by  solving  the  system 

CNwN(t,  q )  -j-  AN{q)wN(t,  q)  =  FN{t,  q) 

(3.3a) 

u;"(0)  =  w% 

(3.36) 

yN{t ,9)  =  J2 

i<N 

(3.4) 

where 

[CN)ij  :=<^.<  >  for  *,J  =  1,2,...,IV 

:=  o{q){<t>" ,$?)  for  i,  j  =  1,2,  •  •  •  ,7V 

[«;*(*,  g)],  :=  w?(t,q)  for  t  =  1,2, •••,7V 

[FN(<,g)]<:=L(<,?)(^)  for  i=l,2, ...,7V 

12 

\  ^  “*«.  **  *m  ^ **■,  m  **k  ^  m  ■**  *  m  ■”*  *  .  W  »  *  .  ■* 


K^],:=<trooT-1(<7),^> 


for  i  =  1, 2,  •  •  • ,  N. 


The  approximating  identification  problems  thus  take  the  following  form: 

(AIDP)n  Find  qNeQ  which  minimizes 

=  (3-5) 

2  Jo 

subject  to  the  approximating  system  (3.3)  and  (3.4). 

Our  convergence  results  for  the  finite  element  schemes  are  summarized  in  the  following 
two  theorems. 

Theorem  2:  Let  C  Q  be  a  sequence  such  that  qM  — >  qeQ  as  M  — *  oo  and 

let  uN(qM)  and  u(q)  be  the  solutions  of  Eqs.  (3.3)  and  (2.7)  corresponding  to  qM  and  q, 
respectively.  Then,  under  hypotheses  (H-0)  to  (H-6),  we  have  uN(qM)  — ►  u{q)  strongly  in 
L2(T ;  Hl(G))  as  N,M-+oo. 

Theorem  3:  Let  qN  be  a  solution  of  the  problem  (AIDP)N .  Then  the  sequence  {gw} 
admits  a  convergent  subsequence  {qNk}  with  qNk  —>  q  as  k  —>  oo.  Moreover,  q  is  a  solution 
of  the  problem  (IDP). 

The  proof  of  Theorem  2  follows  from  the  general  convergence  framework  for  parameter 
identification  problems  given  in  [3]  and  [4].  To  ensure  the  desired  convergence,  it  suffices  to 
show  that  the  sesquilinear  form  o(g)(-,-)  satisfies  the  continuity,  coercivity  and  bounded¬ 
ness  conditions  as  stated  in  [3],  [4],  But  this  is  a  result  of  Theorem  1  under  the  hypotheses 
(H-0)  to  (H-4). 

The  proof  of  Theorem  3  can  be  carried  out  by  using  Theorem  2  and  the  compactness 
of  Q.  Since  qN  is  a  solution  of  the  problem  [AIDP)N ,  it  is  clear  that 

JN{qN)  <  JN{q)  for  MqeQ.  (3.6) 

Thus,  if  we  can  argue  that  for  any  qM  — *  q  in  Q, 

yN{qM)-^y{<l)  in  L2(T;  F)  as  N,M o o, 
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then,  we  can  obtain  the  desired  inequality 


J{q)  <  J{q)  for  VyeQ 


by  taking  limits  in  (3.6).  But  the  needed  arguments  follow  immediately  from  Theorem  2 


\yN(qNk)  y(g)||£2(T;F)  <  K\\uN(qNk)  —  u(q) ||x,a(r;v) 


where  K  is  independent  of  qNk  and  q. 


IV.  OPTIMIZATION  TECHNIQUES  FOR  THE  APPROXIMATE  ESTIMA¬ 


TION  PROBLEMS 


Let  qN  be  an  optimal  solution  of  the  problem  (AIDP)N .  Then  a  necessary  condition 


for  qN  to  be  optimal  is  characterized  by 


VqJN{qN)  .(q-qN)>  0  for  VqeQ 


where  Vq  denotes  the  gradient  of  JN{q)  with  respect  to  q.  From  Eq.  (3.5),  we  have  for 


k  =  1,2,-**  ,n 


=  /0',K(i.«))'(c,V((,9)  -y/W)* 


where 


[C?]i, i  =<  4>?,  4>"  >f  for  i,j  =  1, 2,  •  •  • ,  N 


[YdN(t)}j  =  <  yd(t),4>?  >f  for  j  =  l,2,-..,iV. 


Using  the  same  procedure  as  in  [9],  we  can  evaluate  the  gradient  vector  by  (for  k  = 

1,2,  •••,«) 


where  vN{t,q)  is  the  solution  of  the  adjoint  equation, 


-CNvN(t,  q)  +  A'N(q)v"(t,q)  =  YdN (t)  -  CbNwN(t,q) 


(4.3  a) 


$ 

fy 


vN{tf,q)  =  0. 


(4.36) 


In  Eq.  (4.3a),  the  matrix  A*N(q)  is  given  by 


[^w(«)W  =  *•(«)(*?.*?) 

where  a* (<?)(•,•)  is  the  adjoint  sesquilinear  form  of  a(q)(-,-)  defined  by 

°*  (q)(^>,^) 

:=  /  fd  £  a<i(9»z)^-|7-EM9.*)^|^  +  («0-E|7(9*2))^ 

^  G  L»,J<2  UZjUZt  y<2 

f  *  £/i  f  * 

+  /  “7 - 

Jo  r^,2j)  Jo 

/-1  ^r'^.Zi)  ^r'^Zi)  ,1fJ  ,, 

+/„ 

f1  Cir'{q,Zi)^ 

-  Jo 

Consequently,  the  optimality  condition  (4.1)  of  the  problem  ( AIDP)N  can  be  characterized 


£  I''  vN(t,  $")'{[  qN)  -  V1kFN(t,  qN)}(qk  -  q")  >  0  for  all  qeQ. 

k= l  Jo 

(4.4) 

In  the  sequel,  we  discuss  computer  implementation  of  numerical  schemes  for  the  prob¬ 
lem  (AIDP)n .  Since  we  can  evaluate  the  gradient  of  the  cost  function  using  (4.2),  many 
optimization  techniques  for  the  constrained  problems  are  readily  applicable  to  our  problem 
(see  [11]  and  the  references  therein).  For  ease  in  exposition,  here  the  compact  set  Q  C  RN 


is  assumed  to  be  defined  by 


Q  =  {q  =  {qi,qi,---,qn)eRN\Tiq  <  <?} 


where  II  and  q  denote  a  given  constraint  matrix  and  vector,  respectively.  For  the  numer¬ 
ical  results  reported  in  this  paper,  we  used  the  gradient  projection  method  [12]  which  is 


a  particularly  useful  technique  for  optimization  problems  with  the  linear  inequality  con¬ 
straints  such  as  those  given  in  (4.5).  We  use  this  method  as  presented  in  [12];  the  iterative 
algorithm  for  finding  qN  can  thus  be  stated  as  follows: 


Step  0:  Choose  an  initial  value  q^  in  Q  and  set  i  =  0. 

Step  1:  If  II?W  <  ^  Set 

and  proceed  to  Step  8;  otherwise,  proceed  to  Step  8. 

Step  2:  Compute  the  current  direction  by 

M  _  _PVJJ%I^L 

where 

p  =  /-n;(n„n;)-1np 

and  IIP  includes  the  gradient  of  all  currently  active  constraints  associated  with  matrix  II. 
If  gW  ^  0,  proceed  to  Step  8;  otherwise,  proceed  to  Step  4. 

Step  3:  Compute  A^m  satisfying 

JN(q{i)  +  A J2n*W)  =  min  JN(q^  +  A 

A«(0.AJ 

where  A  is  the  largest  step  that  may  be  taken  from  along  g ^  without  violating  any 
constraint.  If  A^tfl  =  A,  then  add  the  new  contraints  to  the  matrix  ITP  and  proceed  to  Step 
4 ;  otherwise,  the  new  approximation  to  the  solution  is  given  by 

9<i+1)  =  +  A<Lj(0- 

Replace  i  -|-  1  by  t  and  return  to  Step  1. 

Step  4-'  Compute  the  vector  6(q)  by 


*(«w)  =  -(n,n;)-‘nrv,,j(?('>). 


i 

n 


If  all  components  of  0  are  nonnegative,  then  set 


5"  = 


and  terminate  the  computation;  otherwise,  delete  the  column  of  Up  corresponding  to  the 
smallest  component  of  0(q W),  replace  i  +  1  by  i  and  return  to  Step  1. 


V.  NUMERICAL  PROCEDURES 


In  a  series  of  numerical  experiments,  we  used  a  test  example  constructed  as  follows: 
We  chose  a  function  r(q),  generated  the  corresponding  solution  numerically,  added  random 
noise,  and  then  used  this  as  “data”  for  our  inverse  algorithm.  The  parameter  function 
r(£, q)  to  be  identified  is  a  piecewise  cubic  polynomial  function  (see  [6]  for  more  details). 
We  denote  the  knot  sequence  for  r  by 

0  =  r0"  <  t?  <  ■  ■  •  <  r”  <  r"+1  =  1 

and  the  unknown  function  r(£,q)  is  given  by 

r(£.0  =  P?tt) 

=  a-i.i  +  -  t?)  +  a3>,(£  -  t?)2/2  +  a4,,(£  -  t?)3/6  (5.1) 

for  t?  <  £  <  r,"  j  *  =  0,1,  •••,»*. 

The  unknown  parameter  vector  q  =  {<7i}"=i  is  then  given  by 


9,  =  r(n)  for  *  =  1,2,  •••,«. 


Further,  we  assume 


Po(0,?)  =  pn{l,q)  =  t  (5.3) 

Po(°<q)  =  P'M’V)  =  (5.4) 

Substituting  (5.2),  (5.3),  and  (5.4)  into  (5.1),  the  coefficients  {a*,}  can  be  determined 
uniquely  and  r((,q)  satisfies  the  hypotheses  (H-l),  (H-2),  and  (H-4).  In  order  to  guarantee 
the  hypothesis  (H-3),  we  impose  the  constraints 

0i  <  9.  <  02  for  i  =  1,2,  ■■■  ,n. 


'  n+m"  •  ■ 


! 

I 


Hence,  the  matrix  II  (2 n  x  n)  and  the  vector  q  (2 n  x  1)  defining  the  admissible  parameter 
class  Q  (see  (4.5))  is  given  by 


n  = 


i 

-l 


02 

-01 


02 

-01 


To  discretize  the  system  model  by  the  finite  element  method,  the  domain  G  is  divided 
into  a  finite  number  of  elements  {efc}*Li  (K  <  N)  and  a  number  of  nodes  defined  by 
{z,  =  (zj,  ^2 ) } i  are  selected  in  G.  For  convenience  of  computations,  we  set  l  =  1  in  G. 
Each  element  is  preassigned  as  an  axiparallel  rectangle  with  nodes  at  the  vertices.  The 
restriction  of  4>0  to  any  element  e*  is  given  by  the  bilinear  polynomial  form, 


4fc + ciiz  i  +  ci3lz  2  + 


(5.5) 


for  z=(zi,z2)eek  k~l,2,--,K  and  *  =  1, 2,  •  •  • ,  N. 

The  coefficients  {c,-^}  can  be  chosen  such  that  each  polynomial  form  (5.5)  satisfies  the 
properties  of  a  piecewise  bilinear  basis  function  (see  e.g.,  [1],  Ch.  5).  The  integration  of 
element  matrices  CN,  AN(q),  A*N(q),  and  Cf,  and  the  element  vectors  FN  and  Yf  can 
be  computed  numerically  by  a  Gauss-Legendre  formula.  Thus,  the  state  model  (3.3)  and 
its  adjoint  system  (4.3)  can  be  solved  numerically  by  an  implicit  scheme  with  respect  to 
discrete  time  t  =  ih  (t  =  0, 1,  •  •  ■ ,  m)  where  h  =  tf/m.  The  evaluation  of  cost  functional  JN 
and  its  gradient  V,  JN  is  the  computationally  expensive  part  of  our  algorithm  since  these 
involve  the  integration  of  the  states  wN(t,q)  and  the  adjoint  states  vN(t,q )  with  respect 
to  time  t  over  T.  This  can  be  accomplished  by  using  the  two-point  Gauss  formula. 


The  input  data  are  preassigned  as 


The  known  parameters  ct ,  c0,  and  h  in  Eq.  (1.1)  were  set  as 


C\  =  0.034,  c2  =  0.001,  h  =  0.1. 

The  observed  data  {^(f)}  were  generated  by  solving  the  finite  element  model  (3.3).  The 
number  of  finite  elements  and  nodes  in  the  numerical  experiments  were  set  as  K  =  256  (= 
16  x  16)  and  N  =  289  (=  17  x  17),  respectively.  The  final  time  and  number  of  time 
divisions  were  taken  as  tj  =  10  and  m  =  100.  Random  noise  at  various  levels  from  0%  to 
50%  was  added  to  the  numerical  solution,  thereby  producing  simulated  noisy  “data”  for 
the  algorithm.  The  set  £  relative  to  data  acquisition  was  given  by 

£  =  U  N„ 

1=1 

where  NZ]  denotes  a  neighborhood  of  points  x,  at  dG^,  i.e., 

Nx.  —  (i j  —  e,  x j  +  e)  for  j  =  1, 2,  •  •  • ,  p. 

Using  such  data,  the  estimation  algorithm  given  in  Section  4  was  tested. 

Example  1:  In  this  example,  the  dimension  of  unknown  vector  was  taken  as  n  =  4  and  the 
knot  sequence  was  given  by 

t*  =  i/ 5  for  i  =  1, 2,  •  •  • ,  5. 

The  values  of  the  true  parameters  were  chosen  as 

qi  =  r(r,4)  =  0.8  for  i  =  1,2, 3,  4. 

The  lower  and  upper  bounds  of  the  unknown  parameter  vector  were  taken  as  /?j  =  0.3  and 
02  =  1-1,  respectively.  The  initial  guesses  for  the  parameters  were  given  by 

<?J0)  =  1  for  *  =  1,2, 3,4. 

The  number  of  sensors  was  taken  as  p  =  9.  Table  1  shows  the  estimated  parameter  numer¬ 
ical  results  for  the  data  with  noise  free,  5%,  10%,  and  50%  relative  noise  and  Figure  5.1 

19 


shows  the  estimated  parameter  function  r(£,  qN)  and  true  function  r(£,  q)  which  correspond 
to  the  estimated  boundary  shape  and  true  boundary  for  the  10%  noise  case. 


9i 

<72 

h 

<74 

BMWMI 

True  Value 

0.800 

0.800 

0.800 

0.800 

Initial  Guess 

1.000 

1.000 

1.000 

1.000 

Noise 

iteration  6 

0.850 

0.882 

0.880 

0.849 

3.36  x  10~2 

Free 

iteration  13 

0.820 

0.819 

0.821 

0.819 

9.92  x  10"3 

iteration  17 

0.810 

0.785 

0.810 

0.806 

5.39  x  10~3 

5% 

iteration  6 

0.851 

0.879 

0.878 

0.844 

BrprFMM 

Noise 

iteration  13 

0.828 

0.821 

0.818 

0.820 

iteration  17 

0.815 

0.797 

0.791 

0.814 

10% 

iteration  6 

0.849 

0.853 

0.861 

0.834 

Noise 

iteration  14 

0.787 

0.847 

0.835 

0.750 

iteration  18 

0.827 

0.792 

0.799 

0.796 

50% 

iteration  7 

0.922 

0.820 

0.748 

0.808 

3.36  x  10~2 

Noise 

iteration  15 

0.902 

0.793 

0.772 

0.886 

3.40  x  10~2 

iteration  19 

0.813 

0.783 

0.727 

0.794 

1.90  X  10~2 

Table  5.1.  True  Value  and  Estimated  Values  in  Example  1. 


Example  2:  We  chose  the  same  dimension  of  unknown  parameter  vector  as  in  Example  1 
and  we  also  used  the  same  knot  sequence.  In  this  example,  however,  the  values  of  the  true 
parameters  were  preassigned  as 

<h  =  =  0.9 


and 


92  =  93  =  0.6, 


respectively.  The  lower  and  upper  bounds,  initial  guess  of  unknown  vector,  and  number 
of  sensors  were  given  by  the  same  values  as  in  Example  1.  Table  5.2  shows  the  numerical 
results  obtained  here  for  the  various  sets  of  noisy  data.  Figures  5.2  and  5.3  represent  the 
estimated  parameter  function  for  the  case  of  20%  and  50%  noisy  observation  case. 


w 

5»» 


True  Value 
Initial  Guess 
Noise  iteration  16 
Free  iteration  25 
5%  iteration  16 
Noise  iteration  26 
10%  iteration  16 


Noise  iteration  27 


iteration  16 


Noise  iteration  27 


25%  iteration  16 


Noise  iteration  26 


iteration  16 


9i 

92 

93 

94 

0.900 

0.600 

0.600 

0.900 

1.000 

1.000 

1.000 

1.000 

0.817 

0.792 

0.778 

0.821 

0.894 

0.607 

0.602 

0.893 

0.953 

0.694 

0.811 

0.906 

0.896 

0.604 

0.605 

0.898 

0.902 

0.807 

0.674 

0.949 

0.908 

0.594 

0.581 

0.896 

0.917 

0.683 

0.812 

0.904 

0.902 

0.603 

0.610 

0.887 

0.939 

0.819 

0.684 

0.962 

0.868 

0.563 

0.563 

0.874 

1.02 

0.618 

0.680 

0.915 

0.951 

0.574 

0.599 

0.953 

7.14  x 
2.85  x 

5.93  x 

1.93  x 
5.62  x 
5.39  x 
5.70  x 
4.19  x 

6.15  X 
1.66  x 
3.64  X 
1.95  x 


Noise  iteration  28  0.951  0.574  0.599  0.953  1.95  x  10“ 

Table  5.2.  True  Value  and  Estimated  Values  in  Example  2. 


Example  3:  In  this  example,  we  deal  with  a  somewhat  more  difficult  case  as  compared 
with  Examples  1  and  2.  We  set  the  dimension  of  parameter  space  as  n  =  8  and  we  chose 
the  knot  sequence  as 


Mlo  *i=i/9  fort' =  0,1,  2,  •••,9. 


True  parameter  values  were  given  by 


9i  =  9s  =  0.99, 


q2  =  q7  =  0.98, 
93  =  94  =  0.94, 


95  =  96  -  0.60, 


respectively.  Figure  5.4  shows  the  corresponding  boundary  shape  to  be  identified.  The 
number  of  sensors  was  taken  as  p  =  17.  The  bounds  and  initial  guesses  for  the  parameter 


Figure  5.2.  True  Function  and  Estimated  Function  in  Example  2  (20%  Noise) 


Figure  5.3.  True  Function  and  Estimated  Function  in  Example  2  (50%  Noise) 


Figure  5.4.  Unknown  boundary  shape  in  Example  3. 


vector  were  the  same  as  in  Examples  1  and  2.  We  ran  numerical  experiments  for  the  case 
of  noise  free,  5%,  10%,  20%,  and  50%  noisy  observations.  Table  5.3  shows  the  estimated 
parameter  vector  obtained  here.  Figure  5.5  represents  the  estimated  boundary  curve  for 
the  10%  noisy  data. 


m 


Table  5.3.  True  Value  and  Estimated  Values  in  Example  3. 


Throughout  the  numerical  experiments,  we  checked  the  robustness  of  the  algorithm  with 
respect  to  noise  in  the  observed  data.  Results  in  three  examples  indicated  that  the  algo¬ 
rithm  worked  very  well  (i.e.,  as  expected)  for  various  noise  levels.  Furthermore,  we  checked 
the  sensitivity  of  the  algorithm  with  respect  to  the  number  of  sensors.  Specifically,  we  com¬ 
pared  in  Examples  2  and  3  the  number  of  sensors  (p)  with  the  dimension  of  parameter 
space  (n).  In  Example  2,  for  data  with  p  =  5(>  n  =  4),  the  algorithm  still  yields  an  almost 
identical  fit  (to  that  for  p  =  9)  even  in  50%  noise  case  while  the  fit  could  not  be  achieved 
under  the  reduced  observation  case  p  -  3(<  n).  Also,  in  Example  3,  (where  n  =  8)  the 
fit  could  not  be  obtained  with  p  =  3  or  p  =  5,  while  the  algorithm  performed  well  with 
P  =  9  (>  n).  Carrying  out  a  large  number  of  other  numerical  tests  in  addition  to  those 
reported  for  Examples  2  and  3,  we  suggest  that  the  algorithm  requires  a  number  of  sensors 
which  is  at  least  equal  to  the  number  of  dimensions  of  parameter  space,  i.e.,  p  >  n. 


VI.  CONCLUDING  REMARKS 


In  this  paper,  we  have  discussed  techniques  for  estimating  the  system  boundary  shape 
in  two  dimensional  parabolic  systems.  By  using  a  simple  coordinate  transformation  tech¬ 
nique,  the  parabolic  PDE  defined  on  unknown  spatially  varying  domain  was  converted 
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Vi 

fa 

fa 

fa 

fa 

fa 

fa 

fa 

sEt<«  I V.-  -  fa2 

True  Value 

0.990 

0.980 

0.940 

0.600 

0.600 

0.940 

0.980 

0.990 

Initial  Guess 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

Noise 

iteration  16 

1.039 

1.027 

0.887 

0.808 

0.674 

0.926 

1.029 

1.033 

3.08  x  10“2 

Free 

iteration  23 

1.014 

1.007 

0.943 

0.624 

0.615 

0.943 

1.003 

1.017 

7.27  x  10"3 

5% 

iteration  16 

1.042 

1.031 

0.879 

0.749 

0.666 

0.920 

1.049 

1.031 

2.85  x  10^2 

Noise 

iteration  23 

1.023 

0.987 

0.948 

0.626 

0.628 

0.937 

1.002 

1.020 

7.90  x  10“3 

10% 

iteration  16 

1.037 

1.058 

0.899 

0.733 

0.721 

0.881 

1.044 

1.025 

2.82  x  10~2 

Noise 

iteration  24 

1.010 

1.009 

0.955 

0.586 

0.578 

0.948 

0.994 

1.000 

6.20  x  10~3 

20% 

iteration  16 

1.034 

0.989 

0.950 

0.592 

0.595 

0.994 

1.018 

1.061 

1.61  x  10~2 

Noise 

iteration  24 

1.022 

0.985 

0.915 

0.616 

0.610 

0.941 

0.993 

1.061 

1.06  x  10~2 

50% 

iteration  16 

1.090 

1.176 

0.872 

0.637 

0.654 

0.933 

1.166 

1.096 

4.01  x  10“2 

Noise 

iteration  24 

1.030 

1.033 

0.939 

0.566 

0.603 

0.933 

1.046 

1.052 

1.47  x  10~2 

■Pi 
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Figure  5.5.  True  Function  and  Estimated  Function  in  Example  3  (10%  Noise) 
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into  the  same  type  PDE  with  unknown  coefficients  defined  on  a  fixed  domain.  Thus,  our 
fundamental  approach  was  placed  within  the  theoretical  framework  for  parameter  identifi¬ 
cation  problems  given  in  [2], [3],  and  [4].  The  practical  utility  of  our  algorithm  is  supported 
through  a  series  of  numerical  experiments,  a  summary  of  which  is  given  in  Section  5.  These 
simulations  were  carried  out  on  the  Sun  Microsystems  at  ICASE,  NASA  Langley  Research 
Center.  For  three  different  numerical  examples,  using  data  with  no  noise,  the  proposed 
algorithm  yields  an  almost  perfect  fit,  while,  as  expected,  the  fit  degenerates  significantly 
as  noise  in  the  observation  becomes  more  pronounced. 

Although  here  we  discuss  only  the  case  where  the  unknown  boundary  shape  is  rep¬ 
resented  by  a  simple  function  of  one  variable,  our  basic  parameter  estimation  ideas  and 
techniques  can  be  readily  extended  to  consider  more  general  classes  of  geometrical  struc¬ 
tures  for  the  system  boundary.  For  example,  we  may  also  treat  the  case  where  the  unknown 
boundary  shape  is  characterized  by 


r(g,xi,x2)=0  for  (x!,X2)eR2. 


We  are  currently  pursuing  investigations  for  these  cases. 
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